DMRG-study of current and activity fluctuations near non-equilibrium phase 
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Cumulants of a fluctuating current can be obtained from a free energy-like generating function 
which for Markov processes equals the largest eigenvalue of a generalized generator. We determine 
this eigenvalue with the DMRG for stochastic systems. We calculate the variance of the current in 
the different phases, and at the phase transitions, of the totally asymmetric exclusion process. Our 
results can be described in the terms of a scaling ansatz that involves the dynamical exponent z. 
We also calculate the generating function of the activity near the absorbing state transition of the 
contact process. Its scaling properties can be expressed in terms of known critical exponents. 
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PACS numbers: 

Physical systems that are in contact with two reser- 
voirs at a different temperature or chemical potential, 
develop a heat or particle current [1]. In macroscopic 
systems, fluctuations of these currents can often be ne- 
glected. As is the case in equilibrium systems, one can 
however expect that such fluctuations become impor- 
tant in mesoscopic systems and in the vicinity of a non- 
equilibrium critical point • 

The statistics of current fluctuations in mesoscopic 
conductors have received a lot of attention in the past 
decade Q, since they can, for example, give insight on 
correlated electron transport. It is nowadays possible 
to measure experimentally third and higher order cumu- 
lants of the current in problems of charge transport 0, 0| ■ 
Theoretically, these cumulants can be obtained as deriva- 
tives of a generating function. This function has many 
similarities to the free energy in equilibrium systems. 

In the present Letter, we focus on the scaling of the cur- 
rent distribution in one-dimensional (classical) stochastic 
models such as the (a)symmetric exclusion process. This 
stochastic process is a standard model of non-equilibrium 
statistical mechanics 0] ■ Rigorous results are known 
for the current distribution in this model both on a ring 
and for open boundaries [1, Sll^, El- Moreover, several 
approximate and numerical approaches to this problem 
have been developpcd: simulation techniques that sample 
rare events [13, [l3| , renormalisation approaches [ij] and 
perturbation techniques fl^. Here we apply for the first 
time the density matrix renormalisation group (DMRG) 
to the investigation of current fluctuations. We illustrate 
the method for the current of the totally asymmetric ex- 
clusion process, but the technique is more general. As an 
example we also present results on the total number of 
changes of configuration (a quantity that has been called 
activity [3]) in the contact process [2]. 

In the totally asymmetric exclusion process (TASEP), 
each site of a one-dimensional lattice of L sites can be 



empty or occupied by at most one particle. The dy- 
namics of the model is a continuous time Markov pro- 
cess in which a particle hops to its right neighbor with 
unit rate provided that site is empty. At the left bound- 
ary particles enter the system with rate a, while at the 
right boundary they leave it with rate /3. Asymptoti- 
cally, the TASEP reaches a non-equilibrium steady state 
(NESS) in which a nonzero current flows through the 
system. It is by now well known that the TASEP has 
three distinct phases in the low density (LD) phase 
(a < 1/2, (3 > a) the average current J (per bond and in 
the thermodynamic limit) equals a{\ — a) while in the 
high density (HD) phase (^ < 1/2, a > /3) it is /3(1 - (3). 
Finally, in the maximal current (MG) phase, J = 1/4. 

Let JL{t) be the total current through all bonds up to 
time t during a realisation of the stochastic process. The 
statistical properties of this current can be obtained from 
its generating function 



^(s,L) = lim ilog(e"'^^(*)) 



(1) 



where the average is taken over the realisations of the 
stochastic process. The distribution of JL{t) at large 
times can be determined from /i(s, L) by a Legendre 
transformation while the average current J{L), its vari- 
ance A(L) and higher cumulants can be found as deriva- 
tives of /x(s, L): 

J{L) ^ hm hj^(t))^^{Q,L) 

t^oo t OS 



A(i) = hm i((j!(t))-(J40)')-§(0,i)(2) 



t^oo t 



Let ai be a spin variable which equals 1 if the site i is 
vacant and —1 if it is occupied. The state of the system 
is then characterized by the probability P(C; t) to be in 
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a given microstate C = {cti, . . . ,0-^}. This probabiHty 
evolves according to the master equation [181] 
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dP{C; t) 
dt 



= HP{C;t) 



(3) 



where H is the generator of the stochastic process. The 
properties of the NESS of the stochastic process can be 
determined from the (right) eigenvector of H with the 
largest eigenvalue [3]. Similarly, it is not difficult to 
show that generating functions such as /i(s, L) can be 
obtained as the largest eigenvalue of a modified generator 
H{s) [6]. For the current JL{t), H{s) equals 



L-l 



His) = ^ - niVi+i] + a [e^s^ - ni] 

i=l 

+ p [e'sl - vl] 



(4) 



Here we have used the "quantum" notation for stochas- 
tic systems [l^. The operators sf and s~ respectively 
destroy and create a particle at site i, while rii and Vi 
count the number of particles and vacancies at that site. 

Formulated this way, determining ^(s, L) is mathe- 
matically similar to finding the ground state energy of a 
quantum spin or fermion chain. One of the most succes- 
ful numerical techniques to study low temperature prop- 
erties of quantum chains is the DMRG [20|, [21 1. More 



recently, this method has been extended to stochastic 
systems [l^, [l^] where the main difference is that in gen- 
eral the generator H is non-Hermitian. Here we apply 
the method for the first time to operators such as H{s) 
which are neither Hermitian nor stochastic. We found 
that with the DMRG it is possible to obtain /i(s, L) nu- 
merically exact for systems up to L = 60 with only mod- 
est computing facilities. Since there are no essential new 
ingredients in the method as such [2J|, we focus here on 
the results. 

Firstly, in order to test the method we have calculated 
/i(s, L) for the symmetric exclusion process (SEP) for 
which this function is known for large L values [25|. In 
the SEP, particles can hop both to the right and left with 
equal rate. At its boundaries, the system is in contact 
with particle baths of density pa and pb- In Fig. 1, we 
show typical results for L^(s, L) for various L- values, to- 
gether with the exact result (full line). As can be seen, 
there is a fast convergence towards the asymptotic re- 
sults. 

Going back to the TASEP, we calculated fi{s, L) in 
the various regions of the phase diagram. The cumulants 
J{L) and A(L) are then determined by numerical differ- 
entiation. As an example, we present in Fig. 2 our results 
for J(L)/(i + l) and A(L) in the MC phase. Also shown 
are the exact results for JYL) per bond obtained from the 
matrix product ansatz [17| . The numerical data coincide 
with the exact ones within the accuracy. The variance of 
the current was so far not determined exactly. We find 
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FIG. 1: Cumulant generating function for the symmetric ex- 
clusion process with pa = l,Pb = 0. Shown are the asymp- 
totic results of [25l | (full line) and DMRG results for different 
system sizes. 
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FIG. 2: (a) Average current per bond from DMRG (squares) 
compared with exact results from (crosses), (b) Variance 
of the total current. Both results are for a = 3/5 and l3 = 2/3. 



that A(i) increases as L"'. The corrections to this power 
law behavior are strong and cannot be neglected for the 
system sizes we studied. In order to get reliable expo- 
nent estimates we have used the BST-algorithm . We 
find that in the MC phase, but also at the transition line 
between the MC and LD (or HD)-phase, a = 1.50(2). In 
the LD (and HD) phase, a changes to 2.01(4). Finally, 
along the coexistence line between HD and LD phases, 
we find a = 2.03(3). These results strongly suggest that 
M ~ limL^oo A(L)/i^ behaves as an order parameter: 
it is zero in the MC-phase and non-zero in both the LD 
and HD phases. 

Given the similarities between the generating function 
and the free energy, it is natural to ask about the scaling 
properties of p near a phase transition. To focus at- 
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tention, we consider the transition line between LD and 
MC phases (a = 1/2, /3 > 1/2). We propose that under a 
rescahng with a factor b the singular part of /i transforms 
as 



/i(s, Aa, L) ~ b-'fi{by^s, by-Aa, L/b) 



(5) 



where Aa = a — 1/2 and z is the dynamical exponent. 
We conjecture that z replaces the dimension d that ap- 
pears in the scaling of the free energy because fJ.{s,L) 
is a quantity per " unit of time" , whereas the free en- 
ergy is per unit of volume. The exponent ys is a new 
exponent associated with current fluctuations, and z/a is 
like a thermal exponent in equilibrium critical phenom- 
ena. From (12) and it follows that J{L) ^ 2.-^+?^' 
and A(L) jj-^^+'^y^ at the transition. From the exact 
results on J{L) and our data on the variance, we find 
z = 1.50(2) and ys — 1.50(2). This value of z agrees with 
that determined by the Bethe-ansatz [J^l, z — 3/2, thus 
providing strong support to the scaling form ([5]) . We con- 
jecture that also ys = 3/2. Finally, j/q can be obtained 
from dJ{L)/da. This derivative can easily be calculated 
from the exact results, and gives ya = 1/2 [i^ . 

To test our scaling ansatz further, we investigate the 
variance of the current which should scale as 

A(i, Aa) - L-''+^y'H{Ly'^Aa) 

with H a scaling function. To match the numerically de- 
termined behavior of the variance in the different phases, 
H{x) should be constant for a; > 0, and linear in x for 
small X < 0. This implies that M goes to zero linearly 
as the LD-MC transition line is approached from below. 
In Fig. 3, we show our data for A(L, Aa)L^~^^= versus 
]^Va (^^Q. <- g') ^j^g scaling is well satisfied and gives ex- 
ponent values close to those determined above [28|. In a 
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FIG. 3: Scaling plot of A(L, Aa)L^-2^= versus L«°Aa at 
s = 0, Aa < 0,/3 = 2/3. 



similar way, we also checked that the scaling of d^ji/ds^ 
as a function of s at the transition line is well satisfied 
(2^ . Finally, also data on the third cumulant of the cur- 
rent can be well described by our scaling proposal [24 1 . 

As a second application of our approach we study 
the one-dimensional contact process (CP). In this model, 
each site of a lattice can be occupied by at most one par- 
ticle. An occupied site becomes empty with rate 1, while 
an empty site becomes occupied with a rate C,X/2. Here 
^ is the number of occupied neighbors. When A < Ac the 
process reaches an absorbing state in which all sites are 
empty. For A > Ac, and in an infinite system, the model 
reaches a NESS with a finite density p of particles. The 
contact process [2^ is a standard model for phase tran- 
sitions out of an absorbing state It is known from 
extensive numerical investigations that its phase transi- 
tion belongs to the universality class of directed perco- 
lation |3(|. The scaling properties of various quantities 
near Ac are well characterized [i'] . Here we are interested 
in the activity Khit) of the model, which gives the total 
number of changes of configuration in a realization of the 
process up to time t. The generating function of K^lt) 
is 

tt{s,L) = lim - log{e'^^^''>) 

t— *oo t 

This function can again be obtained as the largest eigen- 
value of a generator which in this case equals 



E 



(6^5+ - Hi) - ^{ui^i + ni+l){e^s^ ~ v,) 



(6) 



(no = 71^4-1—0). A finite system will always reach the 
absorbing state asymptotically. To avoid this, we allow 
the creation of particles at the boundary sites. Following 
the reasonings made for the TASEP, we expect that near 
the absorbing state transition, tt scales as 



7r(s, AA,i) = b-'TT{by''s,b^/''^AX,L/b) 



(7) 



Here AA = A — Ac. The exponents z = v\\lvy_ — 1.5805 
and vy^ — 1.09684(6) are known numerically 31] while 
yK is a new exponent. 

It is possible to express yK in terms of other, known, 
exponents. From the dynamics of the model one can 
show that {KL{t)) [11] obeys 



d{KL) 
dt 



i=l 1=1 



d{ni) 
dt 



(8) 



In the NESS, the second term in ([8]) approaches zero, 
whereas the first one becomes equal to 2Lp. The scal- 
ing of p is well known and therefore the average activity 
should scale as 



KiL) 



lim 

t — >-oo 



{KL{t)) 



L^^P/^^FiL^/^^AX) (9) 



4 



Here F is a scaling function and /3 = 0.27649 [3l|. Since 
K{L) is also the first derivative of tt we get from (O and 
dH): yK ^1 + - f3)/iyj_ = 2.3284. We have calculated 
7r(s, AA, L) using the DMRG. In Fig. 4, we show our re- 
sults for the variance of the activity as a function of AA 
and L. At criticality, we find that the average activity 
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FIG. 4: Plot of the variance Ak{L) of the activity of the 
contact proces for (top to bottom) L = 44, 38, 32, 26, 20. 

diverges as L-^^^'^^^ while its variance goes as L'^-^^(^\ 
These exponents are close to —z + uk = 0.7479 and 
—z + 2yK = 3.0763 predicted by the scahng ([7]). Other 
evidence of ([7]) can be seen in Fig. 5 where we present 
a scaling plot of K{s,L) = dn/ds as a function of s at 
AA = 0. This quantity should scale as L'^+V'^ GiL^^ s). 
The numerical data again support this prediction [28[. 
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with the DMRG. We proposed a scaling form for these 
generating functions which is supported by all numerical 
data. We believe that this scaling is quite general and 
can be applied to other models as well. 
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FIG. 5: Scaling plot of dTv/ds{s,AX = 0,L)L'-y'< versus 
sL^^ for 22 < L < 50. 



In summary, we determined the generating function of 
the current in the TASEP and of the activity in the CP 
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